AD-A017  374 


BIOCYBERNETI C FACTORS  IN  HUMAN  PERCEPTION  AND 
MEMORY 

David  C.  Lai 
Stanford  University 


Prepared  for: 

National  Aeronautics  and  Space  Administration 
Advanced  Research  Projects  Agency 


September  1975 


DISTRIBUTED  BY: 

khh 

National  Technical  Information  Service 
U.  S.  DEPARTMENT  OF  COMMERCE 


k 


STAN  FORD  UNIVERSITY 


\ i 

CENTER  FOR  SYSTEMS  RESEARCH 


FINAL  REPORT 


Biocy bernetie  Factors  in  Human 
Perception  and  Memory 


Principal  Investigator 
David  C.  Lai 


Information  Systems  Laboratory 


/ 

4*  ■ 


Reptoducod  by 

NATIONAL  TECHNICAL 
INFORMATION  SERVICE 

US  Department  of  Commerce 
Springfield,  VA  22151 


SEL  75-021 


FINAL  REPORT 

for  the  Advanced  Research  Projects  Agency 
of  the  Department  of  Defense 

B IOCYBERNET 1C  FACTORS  IN  HUMAN  PERCEPTION  AND  MEMORY 


PRINCIPAL  INVESTIGATOR: 
Dr.  David  C.  Lai 


Contract  No.  DAHC 15-72-C-0232 
ARPA  Order  No.  2190 
September  1975 
Technical  Report  No.  6741-5 


DDC 

NOV  19  1975 

IEGEuU  £ 


L-j 


Department  of  Electrical  Engineering 
Stanford  Electronics  Laboratories 
Stanford  University  Stanford,  California  94305 


I 


ID  REWORD 

This  is  the  final  report  of  our  research  on  biocybernetic  factors 
in  human  perception  and  memory  which  was  begun  in  1972  at  Stanford 
University  and  was  supported  partially*  or  fully  by  the  Advanced  Research 
Piojects  Agency  of  the  Department  of  Defense  under  contract  DAHC15-72- 
C-0232,  which  is  now  terminating.  This  report  presents  a brief  summary 
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project. 
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Biocybernet ic  factors  in  Human  Perception  and  Memory 

SUMMARY 

The  objective  of  this  research  project  is  to  develop  biocybernet ic 
techniques  for  use  in  the  analysis  and  development  of  skills  required 
for  the  enhancement  ol  concrete  imager  o„  the  "eidetic"  type.  The  scan 
patterns  of  the  eye  during  inspection  01  scenes  are  treated  as  indicators 
of  the  brain's  strategy  lor  the  intake  of  visual  information.  We 
attempt  to  determine  the  features  that  differentiate  visual  scan  patterns 
associated  with  superior  imagery  from  scan  patterns  associated  with 
inferior  imagery,  and  simultaneously,  to  differentiate  the  EEG  features 
correlated  with  superior  imagery  from  those  correlated  with  inferior 
imagery.  Kor  this  purpose,  we  have  designed  a closely-coupled  man-machine 
system  to  generate  image  enhancement  and  to  train  the  individual  to 
exert  greater  voluntary  conti j1  over  his  own  imagery.  The  models  for 
EEG  signals  and  saccadic  eye  movement  in  the  man-machine  system  have 
been  completed.  These  models  are  used  for  monitoring  and  prediction  of 
EEG  signals  and  eye  posit  ions.  At  this  moment,  all  parts  but  the  feedback 
path  of  the  man-machine  system  have  been  implemented.  We  shall  describe 
in  this  report  the  details  of  these  models  and  dismiss  their  usefulness. 

We  have  assumed  that  the  strategy  of  an  individual  with  superior 
visual  memory  is  to  fixate  his  eyes  at  the  optimal  locations  of  the 
visual  scene  at  the  optimal  time  instants  in  relation  to  EEG.  Through 
the  EEG  model,  we  have  concluded  that  the  EEG  signals  play  the  role 
of  timing  mechanisms  for  visual  information  acquisition  and  processing. 

By  using  our  characterization  of  scan  patterns,  we  are  able  to  show  that 
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the  scan  pattern  is  more  consistent  lor  individuals  with  high  scores 
obtained  from  the  Marks'  Visual-Memory  Task  than  those  with  lower  scores. 
These  results  tend  to  confirm  our  assumption  and  help  guide  us  in 
obtaining  a firmer  grip  of  these  optimal  locations  and  time  instants,  "e 
have  also  developed  and  implemented  a technique  for  presenting  any  part 
of  a scene  at  a specified  foveal  location  at  a specified  instant  of 
time.  If  the  specified  locations  and  time  instants  are  optimal,  we 
expect  that  the  image  will  be  enhanced.  Prediction  techniques  have 
been  developed  for  both  saccadic  eye  movements  and  EEG  waveforms  so 
that  lead  times  may  be  provided  for  the  presentation  of  appropriate 
visual  cues  (stimuli).  We  shall  present  and  discuss  many  of  these 
in  this  report. 
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I . INTRODUCTION 


This  project  has  been  aimed  at  the  development  and  use  of  biocyber- 
netic  concepts  and  techniques  for  analyzing  and  developing  skills  that 
are  essential  for  the  enhancement  of  concrete  images  of  the  "eidetic" 
type.  We  have  concentrated  on  the  problem  of  achieving  biocybernetic 
expansion  of  visual  memory  by  using  a closely-coupled  man-machine  system 
which  performs  real-time  monitoring,  analysis,  and  feedback  of  spatial 
and  temporal  cues  that  serve  as  keys  to  human  memory  encoding  and 
recall.  There  is  strong  evidence  that  these  cues  are  heavily  depended 
on  in  memory  encoding  and  retrieval  by  human  nervous  systems.  The  closely- 
coupled  man-machine  system  is  used  for  the  measurement  and  prediction 
of  human  mnemonic  performance  to  determine  the  spatial  and  temporal 
cues  as  well  as  for  the  control  and  enhancement  of  mnemonic  skills. 

Such  a man-machine  system  for  visual  memory  tracking  and  training  has 
been  designed  and  is  depicted  in  Figure  1.  This  system  is  capable  of 
delivering  optimal  sequences  of  sensory  s; imuiation  conditionally 
related  to  eye  position  and  brain  stat j,  and  thus  to  explore  systematically 
their  relation  to  visual  memory  retrieval. 

It  is  seen  from  the  man-mac  line  system  that  we  emphasize  the  real- 
time monitoring  and  prediction  of  central  nervous  activities  through 
the  EEG  signals  and  through  the  tracking  of  eye  movement  and  eye  fixations. 
Employing  these  real-time  techniques,  we  have  attempted  to  determine 
the  temporal  and  spatial  cues  for  human  memory  encoding  and  recall  by 
monitoring  of  brain  states  and  eye  positions.  Again,  through  this  closely- 
coupled  man-machine  system,  we  utilize  *his  information  to  arrange  the 
desired  coincidences  between  various  brain  states,  eye  positions,  and 
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Eye-movement  measurement 


Response 


Figure  1 Block  diagram  of  the 
closely-coupled  man-machine  system 
for  monitoring  and  training  of  visual  memory  skills 


the  delivery  of  visual  stimulation.  The  visual  stimuli  have  been 
presented  on  projection  screen  or  oscilloscope  displays  for  binocular 
or  monocular  viewing.  The  eye  movement  was  measured  with  the  Biometrics' 
Eye  Movement  Monitor,  Type  SG;  it  is  now  measured  with  the  Stanford 
Research  Institute's  (SRI)  Cornsweet  Eye  Tracker.  The  brain  states 
have  been  monitored  through  the  EEG  signals.  It  is  expected  that 
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greater  control  of  image  persistence  and  image  dissipation  can  be 
obtained  by  combining  these  real-time  monitoring  and  prediction 

techniques  with  a feedback  scheme  to  close  the  control  loop  as  depicted 
in  the  block  diagram  of  the  man-machine  system.  In  other  words,  we 
strive  to  use  the  computer  system  to  supplement  and  strengthen  those 
deficiencies  in  human  memory  that  ordinarily  result  in  image 
dissipation  based  on  the  assumption  that  a superior  scan  pattern  of 
visual  inspection  that  results  in  a superior  memory  is  more  consistent 
and  also  less  probable  of  natural  occurrence  than  a visual  inspection 
strategy  that  is  less  consistent.  To  summarize,  we  attempt  to  steer 
the  subject  toward  improved  encoding  and  decoding  strategies  for  memory 
by  using  the  techniques  which  have  been  deveaoped  and  implemented  on 
this  mar-machine  system. 


The  closely-coupled  man-machine  system  utilizes  computer-based 
models  in  the  system  configuration  to  predict  the  kind  of  stimuli  which 
should  be  used  in  order  to  produce  the  desired  future  responses.  The 
EEG  signals  and  the  eye-movement  measurement  are  shown  as  the  responses 
which  are  closely  controlled  since  the  computer-based  models  are  designed 
to  mimic  the  actual  physiological  processes  in  regard  to  their  stimulus- 
response  relations.  To  be  more  specific,  since  the  visual  cues  (stimuli) 
must  be  presented  at  appropriate  locations  in  the  visual  field  at  the 


right  instants  of  time,  monitoring  and  prediction  on  real-time  basis 
of  both  eye  movement  and  EEG  signals  are  essential.  These  models  have 
been  developed  for  the  purpose  of  monitoring  and  prediction.  We  have 
completed  the  work  on  the  computer-based  models  which  will  be  described 


in  detail  in  the  following  sections. 


These  models  have  also  been 
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implemented  on  the  PDP-15.  All  of  the  man-machine  system,  except  the 
interface  electronic  devices,  can  be  realized  on  the  PDP-15  computer. 

At  the  close  of  this  project,  we  have  implemented  all  but  the  feedback 
path  of  the  man-machine  system. 

In  order  to  determine  the  visual  cues  which  serve  as  keys  to 
both  memory  encoding  and  decoding,  we  have  to  characterize  the  scan 
patterns  of  observers  during  the  inspection  of  scenes.  With  this 
quantitative  characterization,  we  will  then  be  able  to  differentiate 
quantitatively  a superior  scan  pattern  of  visual  inspection  that 
results  in  a superior  memory  from  an  inferior  one.  This  quantitative 
description  is  intended  for  use  in  the  determination  of  those  vital 
visual  cues.  We  have  developed  and  implemented  a statistical  method 
for  the  characterization  of  visual  scan  patterns.  Our  experimental 
results  have  established  that  observers  with  good  visual  memory,  as 
scored  by  Marks'  Visual-Memory  Task,  have  a more  consistent  eye  scan 
pattern  in  quantitative  terms  than  those  with  poor  visual  memory.  The 
characterization  of  the  visual  scan  pattern  and  the  related  eye-movement 
measurement  techniques  will  be  discussed  in  Section  IV. 

Last,  but  not  least,  we  will  present  a technique  for  the  presentation 
of  the  visual  cues  (stimuli)  at  appropriate  locations  in  the  visual 
field  at  the  right  moments.  Since  our  ultimate  goal  of  this  project  is 
t .e  visual  memory  enhancement  through  biocybernetic  techniques,  once  the 
visual  cues  are  determined,  they  can  then  be  presented  to  the  observer 
or  trainee  using  our  technique.  In  other  words,  this  computer-implemented 
technique  will  enable  us  to  deliver  visual  stimulation  to  coincide  with 
the  desired  brain  states  and  eye  positions  since  it  controls,  in  real-time 
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the  monocular  field  of  view  of  an  observer  inspecting  outline  drawings 
of  visual  scenes  on  a graphical  CRT  display.  The  details  will  be  discussed 
in  Section  V. 

In  Section  VI,  we  shall  summarize  our  accomplishments  and  make 
concluding  remarks.  Some  of  these  accomplishments  have  been  published 
or  presented  in  international  scientific  conferences.  The  list  of 
publications  relevant  to  this  research  project  is  included  in 
Appendix  B. 

In  summary,  we  have  designed  a closely-coupled  man-machine  system 
for  the  purpose  of  developing  and  use  of  biocybei netic  techniques  for 
analyzing  and  developing  skills  that  are  essential  for  the  enhancement 
of  concrete  images  of  the  "eidetic"  type.  This  man-machine  system  has 
been  implemented  except  the  feedback  path.  In  this  system,  we  have 
developed  models  for  EEG  signals  during  visual  stimulation  and  lor 
eye  movement  during  inspection  of  stationary  scenes,  respectively.  These 
models  are  then  used  for  the  purpose  of  monitoring  and  prediction  of 
various  brain  states  and  eye  positions  so  that  the  desired  coincidences 
between  various  brain  states,  eye  positions,  and  the  delivery  of 
visual  stimulation  can  be  arranged.  This  arrangement  of  the  desired 
coincidences  can  be  carried  out  by  a technique  we  developed  and  implemented. 
The  electronic  interface  systems  in  the  man-machine  system  serve 
a very  important  role  in  data  acquisition  and  monitoring  of  the  brain 
states  and  eye  positions.  These  electronic  interface  systems  consist 
of  EEG  recording  machines  with  Grass  preamplifiers,  eye-movement  measuring 
instrument  VSRI  Eye  Tracker)  and  related  interfacing  hardware  and 
software.  We  shall  cover  some  of  the  unique  features  of  the  interface 


systems  in  Appendix 
which  includes:  (1) 
to  this  project;  (2) 
are  to  be  published; 


A.  This  report  will  conclude  with  Appendix  B 
a list  of  the  names  of  those  who  have  contributed 
a list  of  all  relevant  papers  that  have  been  or 
and  (3)  a list  of  lectures  and  talks  given. 
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II.  ANALYSIS  AND  MODELING  OF  EYE  MOVEMENTS 
DURING  SCENE  INSPECTION 


In  the  closely-coupled  man-machine  system  depicted  in  Figure  1, 
we  have  two  boxes,  each  labeled  as  computer-based  model;  one  is  for  eye 
movement  and  the  other  is  for  EEG  signals.  These  computer-based  models 
are  used  for  monitoring  ->nd  prediction  of  eye  positions  and  various  brain 
states,  respectively.  The  monitoring  of  the  eye  positions  is  necessary 
for  the  determination  of  the  spatial  cues  that  serve  as  keys  to  visual 
memory  encoding  and  retrieval  and  the  prediction  of  eye  position  is 
essential  for  the  presentation  of  these  cues  at  the  desired  locations 
of  the  visual  field.  In  this  section,  we  shall  discuss  the  development 
of  such  a model  for  eye  movements  during  inspection  of  two-dimensional 
scenes. 

1-  An  Algorithm  for  Automatic  Identification  of  Fixations  and  Saccades 

Eye  movements  during  inspection  of  two-dimensional  stationary  scenes 
consist  mainly  of  two  components:  fixations  and  saccades.  The  first  step 

in  our  analysis  of  eye  movements  is  to  separate  the  data  into  these  two 
components.  We  need  to  determine  the  beginning  and  end  points  of  the 
fixations  and  saccades.  This  helps  us  to  determine  the  sequence  of 
fixations  during  a volitional  scan  of  a two-dimensional  scene.  The 
separation  of  saccades  from  fixations  is  useful  for  the  modeling 
to  be  described  later.  In  the  past,  the  delineation  of  fixations 
and  saccades  and  the  scoring  of  the  fixation  sequences  has  been  done 
by  hand.  We  have  developed  an  algorithm  which  does  this  automatically 
and  can  thus  reduce  large  amounts  of  raw  eye-movement  data  to  a sequence 
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of  fixation  points  very  quickly.  With  modifications,  the  algorithm  can 
perform  in  real  time.  In  addition,  it  gives  us  an  objective  method  of 
data  analysis  as  opposed  to  the  earlier  methods. 

The  detection  of  saccade  onset  is  essentially  a problem  of  signal 
detection.  We  can  consider  the  movements  during  fixation  as  noise  and 
the  saccade  as  the  signal.  At  any  instant,  we  compute  the  quantities 
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where  (x^,y^)  are  samPled  values  of  the  horizontal  and  vertical  components 

of  the  raw  data;  (X  . Y„)  is  the  location  of  the  fixation;  and  [u  ) is 

r f 1 x 

a weighting  function  of  length  k . The  algorithm  indicates  the  onset 
of  a saccade  at  t = 0 if  either  R > T or  R > T , where  T is  a 

x y 

preset  threshold  value;  otherwise,  the  weighting  window  is  moved 

forward.  The  threshold  T is  a parameter  of  the  algorithm  which 
depend  on  the  noise  level  in  the  measurement,  the  size  of  the  smallest 
saccade  that  has  to  be  detected,  the  amount  of  delay  (k)  allowed  before 
detection,  etc.  This  approach  is  essentially  equivalent  to  using  a 
matched  filter  since  the  weights  are  determined  using  the  saccade 

model  to  be  described  later.  Based  on  this  model,  we  can  determine  the 
error  probabilities  for  different  size  saccades  and  set  T accordingly. 

We  use  a slightly  different  approach  to  determine  fixations  since 
the  end  of  a saccade  is  usually  much  slower  than  its  onset.  Now  we 
compute 

n = * jj:l*j  - E^il  < Ts  A Ixj  - E^il  < Ts  , 1 £ j s Jfcj 

where  (g  ) is  an  empirically  determined  weighting  function  and  T is 
i s 

a threshold  value  depending  on  the  noise  variance.  In  effect,  we  are 
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computing  the  number  of  points  in  a sample  window  of  length  Si  which 
falls  within  a square  of  size  T around  the  window  weighted  mean 
(e  i E SiYi|  • The  algorithm  moves  to  the  fixation  state  if  n 

exceeds  a preset  value  ng(<  £)  ; otherwise,  the  sample  window  is  moved 
forward.  In  practice,  a simple  averaging  function  is  found  to  be  adequate 

for  (Ki)  • 

We  have  written  a program  to  implement  this  algorithm.  The  program 
reads  raw  eye-data  from  the  magtape  and  uses  the  algorithm  to  determine 
the  fixation  points.  These  are  plotted  and  numbered.  The  program  also 
outputs  the  starting  and  ending  times  and  the  standard  deviations  of 
each  of  the  fixations.  An  example  is  shown  in  Figures  2,  3,  and  4. 

Figure  2 shows  the  raw  eye-movement  data.  Figure  3 shows  the  sequence 
of  fixation  points  obtained  from  this  algorithm,  and  Figure  4 shows  the 
fixation  points  superimposed  on  the  actual  stimulus  used.  It  is  clear 
that  a considerable  data  reduction  has  been  achieved  and  the  result  is 

easier  to  understand  and  analyze.  Other  advantages  of  this  algorithm 
include  easy  adaptation  to  different  noise  levels  and/or  measuring 
instruments  by  changing  some  parameters,  and  a quantitative  estimation 
of  error  probabilities.  In  fact,  the  algorithm  has  been  used  quite 
successfully  during  the  course  of  this  project  on  data  obtained  from  both 
the  Biometrics  unit  and  the  SUI  eye  tracker  by  simply  changing  the 
thresholds. 

2.  Modeling  of  Saccadic  Eye  Movement 

The  next  step  in  the  analysis  of  eye  movements  is  the  investigation 
of  saccadic  movement.  We  need  to  be  able  to  monitor  and  predict  the 
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Figure  2 
a 20-second  sea 


Raw  eye  movements  for 
n . Scale  - 2.5  deg/inch 
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saccade  duration  and  length  in  order  to  present  visual  cues  (stimuli) 
at  the  appropriate  locations  in  the  visual  field  at  the  desired  momer  s. 
Thus,  the  detailed  dynamics  of  the  saccadic  movement  are  of  interest  to 
us.  Hence,  our  approach  is  first  to  develop  an  appropriate  model  which 
embodies  the  dynamic  characteristics  of  the  saccadic  system.  This  model 
is  then  utilized  for  monitoring  and  prediction  of  the  saccadic  movement. 

The  human  saccadic  eye-movement  control  system  has  been  of  great 
interest  to  researchers  in  neurophysiology  and  bioengineering.  The 
general  shapes  ot  the  position,  velocity,  and  acceleration  curves  of  the 
eye  movement  have  been  reported  extensively.  We  have  attempted  to  fit 
the  model  responses  to  empirical  data  by  classical  least-squares 
techniques.  This  enables  us  to  estimate  certain  parameters  of  the  model. 
These  parameters  which  will  be  used  for  predicting  the  final  eye  positions 
have  only  been  determined  by  indirect  means  in  the  reported  models.  We 
have  also  made  important  modifications  to  the  model  so  that  its  output 
will  accurately  characterize  the  observed  responses.  We  shall  describe 
and  discuss  the  model  modifications,  the  parameter  estimates,  the  curve- 
fitting  results,  and  the  use  in  monitoring  and  prediction. 

lor  the  present,  we  are  concerned  with  horizontal  eye  movements 
only.  ior  this  purpose,  an  experiment  for  the  study  of  horizontal 
saccades  has  been  designed  and  used.  A row  of  LKDs  equidistantly 
located,  only  one  of  which  is  lighted  at  any  given  time,  is  used  as 
stimulus.  We  will  refer  to  the  leftmost  position  as  the  rest  position 
and  the  other  six  as  the  target  positions.  The  target  signal  is  supplied 
by  a digital  computer.  The  subject  is  instructed  to  follow  the  light 
spot,  which  repeatedly  jumps  from  the  rest  position  to  one  of  the  target 
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positions  and  back.  The  successive  target  positions  are  chosen  from  a 
table  of  random  permutations  to  avoid  anticipation  by  the  subjects.  An 
interval  of  at  least  1 second  is  provided  between  each  position  change. 

Eye  movements  are  recorded  and  stored  as  described  elsewhere  in  this 
report . 

Using  the  algorithm  described  above,  we  first  mark  the  starting  points 
of  the  saccades.  Then  the  saccade  position  and  velocity  are  plotted 
and  equal  length  saccades  are  averaged  with  each  sample  lined  up  at  the 
start.  As  an  example,  we  show  position  and  velocity  plots  of  a typical 
saccade  in  Figure  5 and  the  averaged  plots  from  25  saccades  of  one 
subject  in  Figure  6.  The  saccades  observed  agree  with  those  reported 
by  other  researchers  in  all  aspects  such  as  duration,  overshoot, 
response  delay,  etc. 

In  order  to  monitor  and  predict  the  saccadic  eye  movements,  we  intend 
to  fit  the  model  output  to  tho  saccade  position  data  by  estimating  the 
model  parameters  from  the  data.  The  endpoint  of  the  saccade  can  then  be 
predicted  from  the  model  output  by  using  these  parameters. 

Let  the  horizontal  position  x.  at  time  t.  after  the  onset  of 

l l 

the  saccade  be 

Xf  = e(ai,a2,...,ak;ti)  + n.  (1) 

where  a ,a  , ...,a  are  the  model  parameters  and  n is  the  noise 

X ^ K X 

term  which  represents  the  errors  introduced  by  the  instrument,  the  data 
collection  errors,  the  inherent  random  nature  of  the  biological  system, 
etc.  We  shall  use  the  least-squares  estimation  to  determine  a , a , ..., 

X « 
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Figure  6 Averaged  plots  from 
25  saccades  of  10°  amplitude 
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They  are  determined  by  minimizing:  the  mean-square  error 


s - Sj*,  - e<Va2,...ak;t.)) 


where  n is  the  number  of  samples  observed;  i.e.,  by  solving  the  set  of 
equations 

—7- = 0 


dS 


da. 


= 0 


(2) 


dS 


6a. 


= 0 . 


In  the  case  when  g is  a linear  function  of  the  parameters  a ,a  a 

1 2’  ’ k 

tq.  (2)  is  reduced  to  a set  of  linear  equations  which  can  be  solved  to 

obtain  the  parameter  estimates.  If  g is  nonlinear,  then  an  iterative 

process  has  to  be  applied  in  order  to  solve  for  a ,a  ,...a 

12  k 

If  we  assume  the  noise  n.  to  be  white  Gaussian,  the  least  squares 

estimates  are  known  to  be  identical  to  the  maximum  likelihood  estimates 

which  are  unbiased,  consistent,  and  asymptotically  efficient.  With  the 

above  assumption,  we  can  also  estimate  the  varian  s of  th  parameters 
2 

and  use  X -test  and  t-test  for  goodness  of  fit.  This  gives  us  a 
quantitative  measure  of  the  goodness  of  our  model. 

The  functional  form  for  g is  determined  by  the  particular  model 
used.  It  undoubtedly  relates  to  the  dynamics  of  the  eye-movement 
system.  We  have  chosen  to  use  Robinson's  model  for  its  simplicity. 

This  model  is  block  diagramed  in  Figure  7. 
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NI-rMLF  Plant 


Figure  7 A model  characterizing 
dynamic  properties  of  saccadic 
eye  movement 


Here  we  have  shown  the  dynamic  part  of  the  model  only;  the  conti oiler 
part  is  left  cut.  The  input  to  the  system  is  a pulse  of  height  h and 
width  i . This  produced  a saccade  of  amplitude  equal  to  the  area  under 
the  pulse;  i.e.,  ht  . NI  is  the  neural  integrator  and  MLF  the  medial 
longitudinal  fasciculus  which  is  a feed  forward  path  in  parallel  with 
the  integrator  and  provides  a lead  network  to  compensate  for  the  plant 
lag.  R represents  the  oculomotoneuron  firing  rate.  For  a pulse  input, 
we  see  that  R is  a pulse  step  as  desired.  The  plant  is  a second- 
order  over-damped  system  with  time  constants  T^  and  7^  • Robinson 
chose  the  MLF  gain  T = T . However,  he  suggested  that  in  order  to 

J A 

produce  the  observed  overshoots  and  undershoots  T^  > T^  and  T^  < T^  , 

respectively,  should  be  used.  For  our  purpose,  we  have  decided  to 

retain  T,  as  a separate  parameter.  The  output  0 is  the  eye  position. 
3 

The  complete  transfer  function  is  then 


V + 1 

s 


1 

(T  s+l)(T  s+1) 

X w 


For  the  assumed  pulse  input,  the  output  0 (t)  becomes 
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t > 1 


(3) 


The  position  and  velocity  curves  given  by  the  above  for  h = .5  deg/ms 
and  t = 20  ms  corresponding  to  a 10  saccade  are  shown  in  Figure  8. 

The  values  of  the  other  parameters  used  are  T = 150  ms  and  T = 7 ms. 

X Ct 

(These  values  are  suggested  by  Kobinson.)  Curves  are  shown  for  T = 190  ms, 

O 

150  ms,  and  110  ms.  Note  the  overshoot  and  undershoot  for  T > T , 

O X 

and  T < T , respectively. 

J X 

Note  that  there  are  four  unknown  parameters,  h » T2  ' anci  * 

For  estimation,  we  have  to  solve  Eq.  (2)  for  k = 4 . However,  from  the 

theoretical  evaluation  of  the  parameter  variances,  we  found  that  the 

estimate  of  T^  will  have  a very  large  variance.  Thus,  the  estimate 

will  not  be  very  meaningful  and  hence,  a fixed  value  of  T = 150  ms 

will  be  used.  As  will  be  seen  later,  this  does  not  affect  the 

results  significantly.  We  have  also  found  that,  for  reasons  of  numerical 

calculation,  it  is  better  to  estimate  T /t  than  T itself.  Thus, 

three  parameters  a = h , a = T , and  a = T/T  were  estimated  by 

1 2 2 3 3 1 

using  the  data  from  t to  tQ  F i , where  t^  is  the  time  of  the  onset 
of  the  saccade  and  t is  the  input  pulse  width  which  is  estimated  from 
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\ ho  saecade  velocity  by  noting  that  the  velocity  reaches  a maximum  at 
t t-  i as  shown  in  Figure  8.  The  estimated  parameters  are  then 
substituted  into  Rq.  (3 ) to  obtain  the  complete  saecade  curve  and  the 
prediction  ol  the  Final  position  is  made  simply  by  extrapolation. 

Typical  results  i i om  the  parameter  estimation  and  prediction 
schemes  ire  shown  in  Figures  9 and  10.  The  model  output  matches  the 
actual  data  quite  closely.  The  goodness  of  fit  has  been  tested  using  the 
X“-test  on  the  sum  of  the  squares  of  the  residuals,  S , and  the  t-test 
on  the  individual  residuals  at  significance  level  P = .05.  Note  that 
the  estimates  of  T ^ are  10.9  ms  and  13.6  ms  as  opposed  to  the  value 
suggested  by  Robinson.  This  result  is  typical  and  the  mean  value  oi 
Tr  is  about  13.0  ms,  although  it  varies  between  subjects.  It  is, 
however,  significantly  different  from  7 ms.  T /T  >1  is  observed  in 
most  cases,  although  the  difference  is  not  significant  as  far  as 
the  final  value  prediction  is  concerned.  The  values  of  h and  i 
depend  on  the  saecade  size  which  implies  that  the  size  is  controlled 
by  both  the  amplitude  and  the  width  of  the  pulse.  The  exact  relationship 
is  being  investigated. 

In  most  responses,  an  overshoot  of  the  final  position  was  observed. 
This  overshoot  is  not  produced  by  the  model,  as  can  be  seen  in  Figures 
9 and  10.  Even  if  we  vary  T^  , the  overshoot  could  not  be  produced. 

This  is  because  fo’  T > T , the  overshoot  decays  with  a time  constant 
T ^ , whereas  the  actual  overshoot  decays  much  faster  (at  about  the  same 
rate  as  the  initial  rise).  In  addition,  we  could  not  fit  the  model 
response  beyond  t to  the  observed  data  within  the  chosen  significance 
level  P = .05.  This  indicates  that  the  model  must  be  modified  to 
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Figure  8 Model  response  for  T = 150  ms, 
^ 2 ~ 7 ms*  T = 20  ms,  h = .5  deg/ms, 

T3  = 190,  150,  110  ms 
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Figure  9 Actual  and  model  simulated 
saccade  position 

= 20  ms,  T = 150  ms,  T2  = 10.9  ms, 
160.5  ms,  t — 16  ms,  h=  .28  deg/ms 
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Figure  10  Actual  and  model  simulated 
saccade  position 

t = 20  ms,  T,  = 150  ms.  T = 13.6  ms, 
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T = 168.2  ms,  t = 24  ms,  h = .39  deg/ms 
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produce  the  right  response. 


The  fact  that  the  overshoot  decays  with  a 


time  constant  of  about  Tt>  and  the  input  to  the  agonist  and  the  antagonist 
muscles  used  by  Clark  and  Stark  suggests  that  the  input  to  the  model 
consists  of  a combination  of  pulses  rather  than  a single  pulse.  A 
study  of  the  velocity  curve  in  figure  6 also  suggests  that  a positive 
pulse  followed  by  a negative  pulse  can  produce  the  desired  output. 

Thus,  our  model  input  is  changed  to  two  pulses  as  shown  below: 


1 1 

hl  pl-l 

1 

^ 2 • 

h 

2 


where  the  saccadc  amplitude  = h^  + . Wore  generally,  we  can  use  p 

pulses  with  saccade  amplitude  equal  to 

P 

£ hi.  . 

i=l  1 1 

In  practice,  we  found  that  no  more  than  three  pulses  were  ever  necessary. 

The  results  of  this  modification  are  seen  in  figures  11  and  12. 
for  all  multi-pulse  input  fits,  Ty  - T^  = 150  ms  was  assumed.  liguic 
11  shows  the  results  of  estimation  from  a single  pulse  input.  figure  12 
shows  the  same  saccade  with  the  model  output  from  a two-rulse  input 
keeping  other  parameters  the  same  as  before.  The  fit  in  the  second 
case  is  clearly  much  better.  figure  12  also  shows  the  tit  of  the 
velocity  curves  from  the  data  and  from  the  model.  Thus,  excellent  fits 
are  obtained  to  both  the  position  and  velocity  curves  with  this  model, 
figure  13  shows  a much  more  dramatic  result  where  a three-pulse  input 
to  the  model  is  used.  The  justification  for  using  such  an  input  can  be 
clearly  seen  from  the  velocity  curve  which  first  rises  to  a maximum 
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Figure  11  Actual  and  model  simulated 
saccade  position 

tQ  = 20  ms,  = 150.0  ms,  T = 13.1  ms, 

= 150.0  ms,  t = 18  ms,  h = .36  deg/ms 
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Figure  12  Actual  and  model  simulated 
saccade  position  and  velocity 
with  two-pulse  input 
tQ  = 20  ms,  Tx  = 150  ms,  = 13.1  ms, 

T3  = 150  ms,  ti  = 18  ms,  = .36  deg/ms, 
T 2 = 27  ms>  h2  = -,°5  deg/ms 


value,  then  approaches  a steady  positive  value,  followed  by  a rapid 
fall  to  a negative  maximum  and  decay  to  zero.  This  identities  four 
distinct  regions,  suggesting  four  different  input  levels  or  three  pulses 
(one  level  being  zero).  In  fact,  for  the  fitting,  the  pulse  widths  are 
identified  by  the  instants  where  the  velocity  signal  changes  character. 
This  would  be  easier  to  identify  on  a more  noise-free  system  where  the 
second  derivative  or  acceleration  can  be  observed.  However,  it  is 
clear  that  the  input  is  more  complex  than  a single  pulse.  At  present, 
we  are  trying  to  investigate  the  relationships  between  the  heights  and 
widths  of  the  various  input  levels  for  different  saccade  lengths  for  the 
purpose  of  prediction.  Similar  responses  are  observered  for  equal 
length  saccades  which  indicates  that  some  identifiable  relations  exist. 

We  are  also  trying  to  extend  the  results  to  vertical  and  oblique  saccades. 

The  above  approach  is  useful  for  prediction  of  the  final  position 
of  the  saccade,  given  a part  of  the  saccade;  i.e.,  we  can  use  it  for 
prediction  only  after  a certain  time  has  elapsed  since  the  onset  of 
the  saccade.  However,  we  also  need  to  know  when  a saccade  will  occur 
in  response  to  a stimulus.  Normally,  there  is  a delay  of  between  150 
ms  to  250  ms  from  the  presentation  of  the  stimulus  to  the  saccadic 
response.  If  we  can  predict  the  actual  delay,  we  can  predict  the 
course  of  the  eye  movements  following  the  presentation  of  the  stimulus 
by  using  this  delay  in  conjunction  with  the  above  model.  We  decided 
to  use  information  in  the  EEG  alpha  rhythm  to  predict  the  onset  of  the 
saccade.  Gaarder  et  al.  have  shown  that  fixation  saccades  are  initiated 
during  a particular  phase  of  the  alpha.  This  is  in  accordance  with 
Wiener's  hypothesis  that  the  alpha  cycle  serves  as  a clock  which 
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provides  timing  for  the  data  gathering  and  processing  portions  of  the 
visual  system.  This  result  was  obtained  from  microsaccades  occurring 
during  fixation.  We  want  to  determine  if  it  applies  for  the  large 
saccadic  movements  as  well.  If  so,  we  can  use  an  EEG-alpha-wave 
predictive  scheme  to  predict  the  phase  of  the  waveform  and  then  use 
this  information  for  estimating  the  onset  of  the  saccade. 

The  same  experiment  as  described  earlier  for  studying  horizontal 
saccades  and  in  addition  to  record  simultaneously  the  EEC  signals  collected 
as  described  in  subsequent  sections  has  been  performed.  The  EEG  was 
filtered  in  the  alpha  range.  Two  methods  were  used  to  determine  the 
relation  between  the  saccade  onset  and  the  alpha  phase.  First,  the 

quadrant  of  alpha  during  which  the  saccade  started  was  determined. 

2 

The  -<  -test  was  used  to  test  the  distribution  of  saccades  in  each  of  the 
iour  quadrants,  which  was  found  to  be  significantly  different  (p  < .05) 
from  a uniform  distribution.  This  test  is  identical  to  the  one  used  by 
Gaarder  et  al.  In  addition,  the  alpha  was  averaged  with  the  saccade 
onsets  lined  up.  Results  indicate  that  there  is  indeed  a preferable 
phase  of  the  alpha  rhythm  for  onset  of  the  saccade.  This  fact  can  be 
exploited  for  the  prediction  of  the  onset  of  saccades  if  a scheme  can 
be  found  for  predicting  the  EEG  alpha  waveforms;  espeically,  their 
phasic  variations.  Such  a method  will  be  described  in  the  next  section. 
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III.  MODELING  OF  EEG  SIGNALS 
DURING  VISUAL  STIMULATION  AND  PREDICTION 
OF  EEG  ALPHA  WAVEFORMS 


In  order  to  determine  the  temporal  cues  that  are  important  in  visual 
memory,  it  is  necessary  to  understand  the  relation  between  the  temporal 
structure  of  visual  perception  and  the  EEG  alpha  cycle.  To  this  end, 
we  have  investigated  the  changes  produced  in  the  EEG  signal  by  simple 
visual  stimuli.  In  particular,  we  are  modeling  and  investigating  the 
various  EEG  entrainment  phenomena  caused  by  periodic  photic  stimuli 
and  the  phasic  changes  in  EEG  due  to  photic  stimuli.  This  will  enhance 
our  knowledge  in  several  areas:  it  enables  us  to  unify  several  seemingly 

different  EEG  phenomena  under  one  principle,  provides  a model  for  the 
phase-dependency  of  stimulus  efficiency,  and  enables  us  to  study  the 
phase-dependencies  of  onset  times  of  saccades.  A number  of  researchers 
have  suggested  that  the  alpha  cycle  is  phase  locked  with  a clock  signal 
which  provides  the  time  base  for  a sampled-data  system  consisting  of 
the  eye,  its  control  mechanism,  and  the  visual  cortex.  The  evidence 
for  this  has  been  accumulated  slowly;  it  is  known  that  visual  perception 
times  and  reaction  times  are  on  the  order  of  one  alpha  cycle  and  dependent 
upon  the  phase  of  the  stimulus.  n addition,  some  recent  work  has 
shown  that  a saccade  is  most  easily  initiated  on  certain  phases  of  the 
alpha  rhythm,  which  coupled  with  the  fact  that  vision  is  blocked 
during  a saccade,  suggests  that  the  alpha  cycle  delineates  visual  data 
sampling  and  processing  periods. 

The  above  discussion  outlines  the  motivation  for  modeling  the  EEG 
signal  during  simple  visual  stimulation;  the  modeling  efforts  to  date 
have  been  based  on  a simple  nonlinear  oscillator  which  is  amenable  to 
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analysis  and  simulation. 


1.  A Nonlinear  Mathematical  Model  for  Entrainment  of  KEG  Signals  by 

Periodic  Photic  Stimulation 

The  model  is  a van  der  .'ol  oscillator  representing  the  behavior  o' 
a subject's  EEG  during  periodic  visual  stimulation  by  a stroboscope  01 
a set  of  LED  goggles.  The  stimuli  are  delivered  ior  20  seconds  at 
some  fixed  frequency  followed  by  a 10-second  stimulus-olf  period.  This 
sequence  is  repeated  for  each  desired  stimulus  frequency. 

The  van  der  Pol  oscillator  can  be  represented  by 

/ 2\  2 2 
x - p^l  - x jx  + u x = ^E(t) 

where  x(t)  denotes  the  EEG  signal;  Q is  the  unstimulated  alpha 
frequency;  l(t)  is  the  external  excitation  (stimulus);  and  p is  the 
nonlinear  coupling  coefficient. 

When  E ( t ) is  a sinusoid,  a first  approximation  to  the  solution 
x ( t ) can  be  made  by  a classical  technique,  such  as  harmonic  balancing. 
The  possible  solutions  can  be  classified  in  the  following  way  when 
E ( t ) = EySin(  • t i Eo) : 

(1)  natural  oscillation  (E^  is  small): 

x(t)  = A sin(u,^t  + <t>)  ; 

(2)  harmonic  entrainment  (w  - ) : 

x ( t ) = A ( t ) sin(u)^t  -+-  0 ( t ) ) ; 

(3)  mth-order  subhannonic  entrainment  (mc.^  = w^): 

x ( t ) - A ( t)  sin  (-£-  t + *<t))  5 

(4)  nth-order  superharmonic  entrainment  = m^) : 

x ( t ) = A ( t)  sinCnw^t.  + f>(t)); 
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(5)  combined  frequency  oscillations  (all  other  cases): 


x(t)  = AQ(t)  sin(wQt  + ^(t))  + A^(t)  sin^t  + ^(t)). 

Some  initial  data  have  been  collected  with  a set  of  specially 
constructed  LED  goggles  which  provide  a sinusoidally  modulated  intensity, 
and  the  results  agree  well  with  the  above  solution  classes.  However, 
since  the  experimental  stimulus  iS  often  a train  of  flashes,  we  have 
developed  an  analysis  technique  for  approximating  the  solution  of  the 
van  der  Pol  equation  when  E(t)  is  a pulse  train. 

We  assume  that  the  unperturbed  oscillation  (E(t)  = 0)  is  correctly 
represented  by 

x ( t ) = aQ  sin  0 (t) 

where  the  angular  displacement  can  be  written  in  terms  of  frequency  and 
phase: 


e (t)  = to  t + <t>(t)  . 

We  wish  to  examine  the  result  of  E(t)  being  a series  of  impulses  of 
strength  q at  a frequency  u . 

Following  a method  due  to  Blacquiere,  we  can  assess  the  perturbation 
caused  by  one  impulse.  At  an  angular  displacement 

Q = U)  t + $ 

n On 

we  apply  the  nth  impulse.  Its  effect  will  be  seen  directly  in  £(t)  , 
causing  a small  step  change 


t -+At 

r 2 

A x = / x(T)di  = qu)Q 

*n 


over  time  At  . However,  x will  be  continuous,  so 
To  a first  approximation,  we  can  expand  Ax  and 
Ax  = Aa  sin  0 + ai»  cos  0 =0 


Ax  = 0 . 

Ax  as  below: 
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A*  = Aau  cos  e - a u>  A<>sin  Q = qujf 
yj  n u u n 0 

where  a a represents  the  change  in  amplitude  and  A<t>  represents  the 
change  in  phase  of  the  system.  Solving  the  above  equations  simultan- 
eously for  Aa  and  , we  obtain 

A a = qu  cos  0 
0 n 

and 

qu 

0 

^0  = _ sin  g 

ao 

In  figure  14,  we  demonstrate  the  amplitude  and  phase  perturbations  as 
a result  of  an  impulse  falling  on  angular  displacement  @ 


A <t>(0  ) — sin  6 

n a n 

Figure  14  Phase  and  u 

amplitude  perturbations 

The  unperturbed  solution  would  have  followed  the  dashed  line,  and  the 
perturbed  solution  follows  the  solid  line. 

In  order  to  deal  with  a series  of  these  impulses,  Aa  will  be 
assumed  to  be  negligible.  This  is  reasonable  since  a more  detailed 
analysis  shows  that  the  amplitude  change  decays  exponentially  with 
time.  Let  T^  = and  T^  = ~~~  be  the  alpha  and  stimulus  periods, 
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respectively.  Then  referring  again  to  Figure  14,  we  see  that  applying 
the  nth  pulse  at  displacement  e,  results  in  a new  displacement 

9"  +A^(0n)  * Si,lcc  P«lses  are  spaced  Tj  seconds  apart,  the 
(n+l)th  pulse  will  arrive  at  displacement 

6n+l  = 0n  + A0(fin)  + Tl  . 

clearly,  for  entrainment  to  occur,  we  must  have  the  congruence  relation 

6 . , = 9 (modulo  T ) 
n+1  n o' 

or 

®n  + A*(V  + Tx  H 6n  (modulo  TQ) 

Hence 

A5(“n)  ^ T1  = k T0  ’ k = °»  !»  2 

k - 0 , we  have  the  unlikely  situation  of  resetting  the  phase  of 

oscillation;  the  required  phase  change  AO^)  is  very  large.  We  will 

not  consider  this  case  here. 

F°r  k = 1 ’ the  s°lution  is  called  harmonic  entrainment,  and  for 
k k 2 , it  is  called  kth-order  superharmonic  entrainment,  where  k 
represents  the  number  of  alpha  periods  elapsed  per  stimulus  impulse.  It 
is  easy  to  show  that  this  is  stable  near 

= 0 , 

mod  T 

0 

as  one  would  expect  from  the  form  of 

In  a similar  fashion,  for  subharmonic  entrainment  to  occur, 

0n4m  6n  (modul°  T0> 
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or 


n+m- 1 

9 + E iM0.)  + mT  = J 

n i In 

i=n 


(modulo  T ) 
0 


Hence, 

n+m-1 

S A*(9.)  + mT1  = Tq  . 

i = n 

th 

When  m = 1 we  have  harmonic  entrainment  again,  and  for  m s 2 , m 
order  subharmonic  entrainment  occurs.  The  stability  of  this  solution 
is  currently  being  investigated,  but  it  appears  that  for  m = 2 a 
small  stable  range  exists. 

In  summary,  when  the  stimulus  frequency  is  near  the  alpha  frequency, 
harmonic  entrainment  is  possible.  When  the  alpha  frequency  is  near  an 
integer  multiple  of  the  stimulus  frequency,  superharmonic  entrainment 
is  possible.  And  when  the  stimulus  frequency  is  near  an  integer  multiple 
of  the  alpha  frequency,  subharmonic  entrainment  is  possible. 

When  the  stimulus  frequency  is  not  in  one  of  the  above  ranges, 
combined  frequency  oscillations  exist.  Thus,  we  have  the  same  five 
solution  classes  for  both  sinusoidal  and  impulsive  stimuli.  The  major 
differences  lie  in  *he  amplitudes  required  to  produce  the  various  effects, 
and  in  the  fact  that  the  impulsive  stimulus  will  produce  integer  harmonics 
of  its  fundamental  frequency  due  to  its  inherent  harmonic  content. 

Results  obtained  from  EEG  data  on  the  model  are  described  and  compared 
in  the  next  section. 

2.  Model  Simulation 

It  is  necessary  to  simulate  the  nonlinear  oscillator  model  to 
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obtain  reasonably  accurate  solutions  since  nonlinear  analysis  is  by 
nature  approximate.  We  are  using  a variable  order  diiferential  equation 
solver  (DVDQ)  which  has  an  Adams-Faulkner  predictor  and  an  Adams- 
Moulton  corrector.  It  provides  lor  up  to  twentieth  order  interpolation 
of  output  values  positioned  independently  of  the  changing  stepsize, 
which  greatly  facilitates  signal  processing  of  the  simulated  solution 
by  allowing  a constant  sample  period. 

The  sine  wave  stimulus  is  handled  in  a straightforward  fashion,  but 

the  impulse  train  is  treated  specially.  Since  it  is  known  analytically 

2 

that  the  exact  change  in  x(t)  is  q (by  integration)  at  the  time 
an  impulse  is  called  for,  the  state  vector  is  perturbed  and  the  integration 
restarted.  The  management  of  impulse  arrival  and  sample  output  is 
performed  by  a discrete  event-queuing  scheduler  on  a priority  basis; 
impulse  arrival  has  priority  over  all  events  to  avoid  frequency  shifts 
in  the  impulse  train.  Error  propagation  as  a result  of  restarting 
the  integration  is  linear  with  the  requested  error  tolerance  and 
adequately  small, 

figures  15  and  16  show  the  phase  and  amplitude  perturbations  caused 
by  a single  flash  landing  on  phases  of  0 , n/2  , n , and  3rt/2  radians. 

The  simulations  match  the  first  order  approximations  developed  in  the 
previous  section  very  well.  These  perturbations  produce  appropriate 
entrainment  phenomena,  as  seen  in  the  spectra  of  Figures  17  and  18; 
there  is  one  spectrum  for  each  stimulus  frequency,  representing  5.12 
seconds  of  impulse  driven  simulation  at  a 100-Hz  sampling  rate.  The 
five  phenomena  predicted  in  the  previous  section  can  be  clearly  seen; 
viz.,  harmonic  entrainment,  subharmonic  entrainment,  superharmonic 
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figure  15  Simulated  EEG  perturbations  resulting 
from  flash  at  phases  of  0 and  n/2  radians 
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Figure  16  Simulated  EFG  perturbations  resulting 
from  flash  at  phases  of  « and  3 «/ 2 radians 
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FREQUENCY  IN  HZ 

Figure  17  Simulated  KEG  spectra  for  flash 
frequencies  from  .5  Hz  to  12.5  Hz 


t\  ts 
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figure  18  Simulated  EEG  spectra  for  flash 
frequencies  from  13  to  25  Hz 


40 


entrainment,  combined  frequency  oscillations,  and  stimulus  harmonics, 
lhese  effects  can  again  be  seen  in  Figures  19  and  20  which  are  spectra 
Horn  b.  12  seconds  of  sine  wave  driven  simulations.  These  spectra  will 
be  compared  with  those  obtained  from  experimental  data  in  a section 
below;  the  next  secti.  i discusses  data  collection  and  processing. 

Data  Collection  and  Processing 

The  data  collected  to  test  the  model  have  been  obtained  from  a 
very  simple  experimental  paradigm.  Subjects  are  presented  a series  of 
•stroboscope  flashes  or  sinusoidally  modulated  intensities  through  closed 
eyelids;  the  stimulus  is  delivered  for  20  seconds  at  some  fixed  frequency, 
followed  by  a 10-second  rest  period.  This  sequence  is  repeated  for  each 
desired  frequency,  typically  .5  Hz  through  25  Hz  in  .5  Hz  steps.  The 
irequency  pattern  can  be  selected  arbitrarily.  The  EEG  signals  are 
obtained  from  left  and  right  occipital  electrodes  referenced  to  yoked 
earlobes;  the  ground  electrode  is  on  the  mastoid.  The  Grass  preamp 
bandpasses  the  signal  between  1 Hz  and  300  Hz,  and  the  A/D  analog 
prefilter  is  a 6-pole  Bessel  with  a 100  Hz  lowpass  cutoff  frequency. 
Digitizing  is  at  1000  Hz,  and  subsequent  filtering  is  done  by  a trans- 
versal filter  to  avoid  further  phase  distortion.  Power  spectra  are 
computed  by  East  Fourier  Transform  on  5.12  seconds  of  the  data,  forming 
the  power  spectrum  and  then  averaging  three  such  spectra  4.0  seconds 
after  the  start  of  the  stimulus.  This  procedure  substantially  reduces 
the  noise,  which  is  further  reduced  for  viewing  aesthetics  by  the  non- 
causal  recursive  filter 

H(z)  = — — ~z  + j_ 

4z 
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Figure  19  Simulated  EEG  spectra  for  sine 
frequencies  from  .5  Hz  to  12.5  Hz 
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figure  20  Simulated  KEG  spectra  for  sine 
frequencies  from  13  Hz  to  25  Hz 
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The  resulting  spectra  may  be  superposed  with  hidden  lines  removed  by 
classical  graphics  technics;  a snail  a„„„„t  of  informatlon  ,ust  be 
saved  after  each  spectrum  is  plotted. 

The  stroboscopic  stinuli  are  produced  by  a Crass  PS-2  strobe  lanp 
in  a sound-proof  enclosure.  It  is  driven  by  a relay  contact.  The 
sinusoidally  nodulated  Intensity  stinuli  are  produced  by  a pair  of 
henispherical  goggles,  each  containing  ten  LEDs,  driven  by  a current 
source  controlled  by  a D/A  converter  on  the  POP-15.  The  voltage  supplied 
IS  proportional  to  cap  (sin  v>  since  the  LED  anplitude  is  linear 
with  current  and  sinusoidal  intensity  Is  desired.  Both  stinuli  are 

delivered  through  closed  evelids  to 

eyelids  to  ensure  a uniform  visual  field. 

4-  prison  of  Model  Simulations  and  Exoer.nu.nt,,,  n,„ 

The  comparison  of  model  simulated  results  with  experimental  results 
is  subject  to  difficulties  arising  from  the  noise  in  the  data  and  the 
random  nature  of  „„ch  of  the  EEC  signal.  The  alpha  frepuency  is  not 
•ixed  and  ,s  in  fact  nonstationary;  it  Is  influenced  by  the  subject's 
level  of  alertness  and  other  uncontrollable  factors.  The  response 
of  the  EEC  to  a stimulus  Is  only  measurable  in  the  average,  since 

individual  responses  are  .asked  or  corrupted  by  noise.  Averaging  can 

also  mask  certain  characterisHr*  ... 

naractenstics  of  the  response,  depending  upon  the 

scheme  used. 

The  best  comparison  of  model-produced  results  and  data  for  entrainment 
Phenomena  is  a comparison  of  their  peer  spectra.  This  allows  the  data 
spectrum  to  be  relatively  free  of  noise  by  averaging  successive  spectra 
fro.  the  same  stimulus  epoch.  Examples  of  stroboscopically  stimulated 
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data  spectra  may  be  seen  in  Figures  21  through  24.  Note  the  excellent 
match  overall  with  the  simulated  spectra  discussed  previously.  The 
last  subject  has  virtually  no  alpha  rhythm,  but  exhibits  all  of  the 
entrainment  phenomena  predicted,  so  the  model  accounts  surprisingly  well 
for  extreme  cases,  Fi  ares  25  and  26  show  spectra  for  a sinusoidally 
stimulated  subject.  Note  the  reduction  in  harmonic  amplitude  as  predicted 
by  the  simulation.  Some  detailed  comparisons  of  each  type  of  phenomena 
follow.  Figure  27  shows  the  unstimulated  case  with  the  model  parameters 
adjusted  to  match  the  subject  (records  were  selected  from  Figures  21 
and  22).  Figure  28  illustrates  harmonic  entrainment,  with  subharmonic 
and  superharmonic  entrainment  treated  in  Figures  29  and  30.  ligure  31 
is  an  example  of  combined  frequency  oscillation.  The  model  simulations 
agree  well  with  the  data.  Additional  comparisons  may  be  found  in 
previous  reports. 

Further  comparisons  are  in  progress,  especially  the  comparison  of 
phase  and  amplitude  changes  induced  by  the  stimuli.  Preliminary  results 
indicate  that  the  model  predictions  are  reasonable  first  order  approxima- 
tions to  the  data,  but  that  revisions  in  the  model  will  likely  be 
necessary  to  fully  account  for  the  observed  phenomena.  The  excellence 
of  match  in  the  frequency  domain  bears  out  the  goodness  of  the  phase 
shift  predictions  for  small  phase  shifts,  but  it  appears  that  larger 
shifts  are  occurring  than  predicted  for  certain  cases. 

In  sum,  a nonlinear  model  for  the  behavior  of  the  EEG  during  visual 
stimulation  has  been  analyzed  and  compared  with  actual  EEG  data.  The 
model  accounts  for  several  phenomena  well  in  a unified  fashion,  and 
suggests  other  phenomena  of  interest.  The  model  specifies  a trigonometric 
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Figure  21  EEG  spectra  from  subject  L.  M.  for 
flash  frequencies  .5  Hz  to  12.5  Hz 
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Figure  23  EEC  spectra  from  subject  i.  M. 
flash  frequencies  .5  Hz  to  12.5  Hz 
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Figure  26  EEG  spectra  from  subject  R, 
sine  frequencies  13  Hz  to  25  Hz 
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Fiuure  27  Actual  and  simulated  EEC  spectra 
in  the  absence  of  stimulation 


e 28  Actual  and  simulated  EEG  spectra 
during  flash  stimulation  of  10.5  Hz 


z««snoAn  ni  u3M0d  1333 


55 


gure  30  Actual  and  simuln 
during  flash  stimulation 


Figure  31  Actual  and  simulated  EEG  spectra 
during  flash  stimulation  of  7.5  Hz 


form  lor  the  excitability  of  the  EEG  as  a function  of  phase,  which 
explains  well  the  various  entrainment  phenomena  seen  in  EEG  data. 

Detailed  phase  analyses  are  in  progress  which  indicate  that  the  model 
prediction  is  good  to  a first  degree,  but  needs  some  revision.  One 
possibility  at  this  point  is  the  existence  of  an  excitability  function 
at  half  the  alpha  frequency.  The  phase-dependency  of  the  model  response 
has  considerable  bearing  on  the  phase-dependency  of  saccade  onset. 

5.  Prediction  of  EEG  Alpha  Waveforms 

In  order  to  present  the  visual  stimuli  to  coincide  with  the  most 
favorable  brain  state  and  to  predict  the  onset  of  saccadic  eye  movement, 
it  is  necessary  to  have  the  capability  to  predict  the  EEG  alpha  waveforms, 
in  particular,  their  phasic  features.  In  conjunction  with  the  EEG  model, 
we  have  developed  and  implemented  an  EEG-waveform  predictive  scheme  by 
using  autoregressive  processes.  The  development  and  implementation  of 
this  scheme  will  be  described  in  detail.  Our  discussion  will  be 
centered  around  the  modeling  ol  EEG  data  by  an  autoregressive  process 
and  its  use  for  forecasting. 

The  j itting  of  time-series  models  to  EEG  signals  has  been  treated 
in  the  literature.  They  have  been  specifically  applied  to  EEG  analysis 
in  the  following  areas:  use  for  EEG  spectral  analysis;  use  of  the 
mixed  model  to  define  certain  parameters  for  describing  the  stationary 
parts  of  the  EEG  signal;  and  use  of  the  one-step  ahead  prediction 
error  to  compare  two  different  models.  Our  interest  lies  in  the  prediction 
of  the  alpha  activity  in  the  EEG  for  a longer  lead  time. 

Although  autoregressive  processes  have  been  used  for  prediction 
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of  time  series  in  general,  they  have  not  been  used  for  prediction  of 
EhG  data.  In  those  general  cases,  the  signal  or  time  series  is  assumed 


to  be  stationary.  EEG  signals  are,  in  general,  nonstationary . 
expedience,  they  may  be  considered  p iecewise-stationary ; i.e., 
over  short  intervals  of  time.  Hence,  in  our  application,  we  are 
restricted  to  limited  sample  lengths.  This  restriction  creates 
problems  which  we  shall  discuss  later. 

To  model  a discrete  time  series  X.  , i = 1,  2,  ...,  by 
regressive  process  of  order  p , we  may  write 


an 


x - a X + a X + ...  + a X + e 
t 1 t-1  2 t-2  p t-p  t 


Kor 

stationary 

certain 

auto- 


where  a^,a^,...,  a^  are  the  parameters  of  the  process  to  be  estimated 
and  e^  is  a zero-mean  white  Gaussian  noise.  It  is  generally  assumed 
that  the  process  X^  has  zero  mean;  otherwise,  the  mean  is  subtracted 
from  the  entire  series.  In  general,  the  autocorrelation  function  of  an 
autoregressive  process  consists  of  two  components;  viz.,  decaying 
exponential  and  damped  sinusoidal  waveforms.  It  has  been  recognized 
that  the  real  roots  of  the  characteristic  equation  of  the  autoregressive 
process,  i.e., 

2 n 

! _ a B - a B - . . . - a Bp  = 0 , 

1 l p 

give  rise  to  the  decaying  expoentials  and  its  complex-conjugate-root 
pairs  are  responsible  for  the  damped  sinusoids  in  the  autocorrelation 
function.  We  show  in  Figure  32  the  autocorrelation  function  of  an  EEG 
sample  filtered  in  the  alpha  range.  It  is  seen  that  the  EEG  autocorrelation 
function  may  be  well  approximated  by  a sum  of  decaying  exponentials  and 
damped  sinusoids.  In  Figure  32  , we  superimpose  the  autocorrelation 
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. autocorrelation 


Figure  32  Autocorrelograms:  solid 

curve  from  a typical  EEG  alpha  waveform  and 
dots  from  an  autoregressive  process  fitted 
to  the  data 
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function  of  an  autoregressive  process  fitted  to  EEG  data  to  show  our 
assertion.  This  leads  naturally  to  the  temptation  of  modeling  EEG  signals 
by  autoregressive  processes.  Since  autocorrelation  function  is  not 
a unique  process,  there  ire  other  processes  which  will  give  rise  to 
the  same  autocorrelation  function.  We  chose  an  autoregressive  process 
based  on  its  simplicity  in  implementation. 

To  accomplish  our  aim  in  prediction,  we  have  broken  down  the  problem 
into  two  parts: 

(1)  Estimation  of  the  autoregressive  parameters; 

(2)  Generation  of  "good'  forecasts  according  to  some  criteria 

using  the  above  estimation. 

Since  the  EEG  process  is  considered  piecewise  stationary,  the  parameter 
estimates  have  to  be  updated  continuously.  Thus,  a scheme  is  needed 
for  updating  the  parameter  estimates  as  new  data  arrive.  Eor  estimation 
of  the  parameters,  we  have  used  the  least-squares  principle.  These 
estimated  parameters  are  then  used  to  generate  the  minimum  mean-square 
error  forecasts  for  a specified  lead  time.  We  developed  a new  modified 
scheme  for  this  purpose. 

To  furnish  the  basis  for  our  modified  scheme,  we  shall  describe  a 
scheme  ordinarily  used.  The  development  and  implementation  of  a modified 
scheme  to  suit  our  needs  will  then  be  presented.  Results  have  been 
obtained  by  applying  our  scheme  to  real  EEG  data.  We  will  present  these 
results  to  show  the  utility  of  our  prediction  method. 

We  will  concentrate  our  discussion  to  a parameter  estimation 
method.  Later  we  will  present  our  modifications  to  this  method  in 
order  to  improve  the  performance  for  our  purpose.  Lot  us  consider 
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the  p -order  process 


X = a X + a X + 
t 1 t-1  2 t-2 


+ a X + e 
P t-p  t 


where  e^  is  white  Gaussian  noise.  The  order  of  the  autoregressive 

process  is  estimated  by  successively  fitting  models  with  higher  order 

to  the  data  and  determining  when  the  last  parameter  a becomes  small 

P 

enough. 

The  least-square  estimates  of  the  parameters  can  be  shown  to  satisfy 
the  following  set  ol  linear  equations 
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where  the  Q 's  are  given  by 
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Q . . = XX  + X X . , , + ...  + X X 
i.)  i j 1+1  j+1  n-jrl  n-i+1 
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Knowing  Q^'s  » we  solve  the  set  of  equations  to  obtain  the  least- 
square  estimates  of  the  parameters  a , a , ...,  a . The  least- 

12  p 

square  estimate  of  the  noise  variance  is  given  by 


rr  — 
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n (QU  + |Ql2,Q13’"*,Ql,p+l|  | V V * * * ’ 3p  | ) • 

Previous  discussions  furnish  us  a numerical  procedure  for  estimation  of 

the  parameters  a , a , ...,  a and  n in  the  least-square  sense  from 
1 2 p ue 

the  data  (X^,  X^,  ...,  X^)  , Based  on  this  procedure,  one  can  model 
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a weakly  stationary  process  by  a pth-order  autoregressive  process. 
Since  EEG  data  exhibit  nonstationarity , this  method  has  to  be  modified 
for  our  problem.  Besides,  our  interest  is  in  the  use  of  the  parameter 
for  prediction  rather  than  modeling. 

Having  obtained  the  estimates  of  the  parameters,  we  turn  our 
attention  to  the  problem  of  prediction.  We  need  a scheme  whereby  the 
prediction  can  be  made  for  a certain  lead  time  and  then  when  a new 
data  point  becomes  available,  the  forecast  can  be  updated  without 
repeating  the  original  process  all  over  again. 

We  will  consider  here  minimum  mean-square  forecasts  only.  Denote 
the  estimation  of  X obtained  at  time  t by  X , . At  time 

I t 

t , the  values  of  the  samples  up  to  time  t are  known.  The  minimum 
mean-square  error  forecast  is  given  by  (using  the  projection  theorem): 

Xt+*|t  = aiXt+£|t  + a2Xt+jfc-2|t  + + apXt+A-p|t  * 

which  provides  the  forecasts  in  the  form  of  a difference  equation. 

Given  the  sample  values  Xt  , X^ X^^  , we  can  progressively 

calculate  t » ^ti2|t  ’ ***  using  the  above  equation.  Using 

this  procedure,  we  can  generate  forecasts  for  any  lead  time  utilizing 
the  parameters  estimated  earlier.  However,  the  forecasts  at  each  point 
have  to  be  updated  as  soon  as  the  sample  Xt+1  becomes  available. 

It  is  rather  cumbersome  to  do  in  the  above  form.  Hence,  a different 
approach  is  required. 

Instead  of  writing  the  original  equation  as 

Xt  = aiXt-l  + Vt-2  + *••  + ap\_p  + et  , 

we  will  express  X in  a weighted  infinite  sum  of  e , e . as 
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xt  = E Yt-J  - 

J=o  J J 


where  the  b^'s  can  expressed  as 


bo  - 1 

b,  = Vo 

b2  * aibl  + Vo 


bj  = “lbj-l  + V.i-2  + "•  + apbj-p  - ^ 2 0 • 

The  minimum-mean-square-error  forecast  is 


(4) 


This  form  lends  itself  easily  to  updating.  This  can  be  demonstrated 
by  Xt+£|  t+1  after  Xt+1  arrives.  We  may  write 


and  hence 


Xt4-X|  t-t-l  j=J_1bjet4i-j 


Xt+ilt+l  Xt+i|t  b^-l6t+l 


(5) 


Thus,  at  time  t + 1 , the  forecasts  for  lead  lengths  up  to  £ - 1 
can  be  easily  obtained  from  the  predicted  values  at  time  t for  lead 
lengths  up  to  £ by  using  Eq,  (5)  and  by  recognizing 

6t+l  “ Xt+1  ~ Xt+l|t  ' 

This  provides  a convenient  scheme  for  updating  forecasts  as  more 
data  become  available. 

The  variance  of  the  forecast  error  may  be  easily  obtained  from 
Eq.  (4).  Since  the  £-step  ahead  forecast  error  at  time  t is  given  b 
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et+il t xt+i  " xt+i!t 


i-1 

= j?o  bJGt+i-J  * 


The  variance  V(i)  may  then  be  expressed  as 

i_1  2 2 

v(i)  = E b“  o 

j=o  J c 


(6) 


This  variance  gives  us  a measure  of  the  goodness  of  fit  for  the  model. 
Since  we  have  assumed  the  e's  to  be  normal,  the  conditional  probability 
density  function  p<xt+ilxt*  xt_1  ...  ) is  normal  with  mean  xt+i | t 
and  variance  V(.£)  . Based  on  this  information,  we  can  obtain  confidence 
intervals  for  the  forecasts.  For  example,  95$  of  the  normal  distribution 
mass  lies  within  + 1.96  a of  the  mean  where  o is  the  standard  deviation. 
Hence,  xt+i;|t  + 1.96  x/v”  is  the  95$  confidence  interval  for  the  forecast, 
i.e.,  the  probability  that  the  actual  value  will  lie  in  that  range  is 
.95. 


The  previous  formulation  gives  the  forecasting  procedure  to  predict 
values  i samples  ahead  from  any  time  t . In  order  to  increase  the 
lead  time  for  forecasting,  one  has  to  increase  the  value  of  £ , It  is, 
however,  evident  from  Eq.  (6)  that  the  error  variance  becomes  larger 
as  £ increases  and  the  confidence  interval  for  a given  level  of 
confidence  will  also  increase  accordingly.  Another  way  to  increase 
the  lead  time  would  be  to  increase  the  sampling  interval  and  thus  keep 
the  value  of  i as  low  as  possible  so  that  the  level  of  performance 
of  the  predictor  may  be  maintained.  However,  the  sampling  rate  or 
sampling  interval  for  a signal  is  determined  by  the  bandwidth  of  the 
signal.  If  the  signal  is  furnished  in  digital  form,  then  one  has 
little  choice  in  altering  the  sampling  interval.  One  could,  of  course, 
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increase  this  interval  by  dropping  every  kth  block  of  data  points.  In 
any  event,  the  increase  in  sampling  interval  means  reduction  in  the 
number  of  data  points  in  a given  interval  of  time.  This  reduction  in 
data  points  leads  to  a higher  variance  for  the  estimates.  Besides, 
the  number  of  data  points  may  be  so  much  reduced,  since  the  EEG  signal 
is  considered  stationary  only  for  a short  period  of  time,  that  the 
estimates  are  not  statistically  significant.  For  these  reasons,  we 
reformulate  the  autoregressive  model  in  a more  general  form  as 
Xt  = aiX.  , + aoX*  ol.  + ...  + a x + e 

* 1 t-k  2 t-2k  p t-pk  t 

which,  when  k = 1 , reduces  to  the  case  discussed  earlier.  This  is 
equivalent  to  an  autoregressive  model  for  the  time  series 

(Xt*  Xt-k’  Xt-2k  **•  } ‘ 

The  least-square  estimates  of  the  parameters  can  be  obtained  by  a 


logical  extension  of 

the  earlier 

results ; 
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aiQ2,p+l  B2Q3,p+l  + + apQp+l,p+l  _Ql,p+l 


where  the  Q^jfs  are  now  redefined  as 
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Qu  - 2 xi 
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Qij  " X(i-l)k+lX(j-l)k+l  + + Xn-(j-l)kXn-(i-l)k  * 

i 4 1 , j i 1 . 

The  least-square  estimates  are  the  solutions  of  this  new  set  of  equations. 
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The  forecast  equations  are  modifield  in  a similar  manner  to  give 
Xt+ki|t  ' aiXt+k(i-i) 1 1 + a2Xt+kU-2)j  t + •** 

+ apXt+k(i-p) | t 

and  Fq.  (d)  becomes 

OC 

Xt+ki|t  " ^VtFkCi-j) 

with  the  b^'s  defined  in  terms  of  the  * s as  before.  In  essence, 
the  scheme  processes  the  available  data  in  such  a regrouped  fashion  by 
resampling  the  data  in  k multiples  of  the  original  sampling  interval 
that  the  level  oi  performance  is  not  degraded.  The  performance  of  this 
scheme  on  KEG  data  will  be  discussed  next.  The  usefulness  of  the  scheme 
is  demonstrated  by  our  results. 

We  developed  a computer  program  to  implement  an  algorithm  for 
recursive  estimation  of  parameters  and  forecast  of  future  values  by 
using  the  new  modified  scheme.  The  program  was  written  in  FORTRAN  for 
a POP-15  computer.  As  an  example  for  illustration,  we  applied  our 
scheme  to  real  EEG  data,  which  were  recorded  while  the  subjects  with 
closed  eyes  were  stimulated  by  stroboscopic  flashes  for  50  seconds  at 
the  rate  oi  10  flashes  per  second,  then  no  stimulation  for  50  seconds, 
then  another  50  seconds  of  stimulation,  etc.  The  signals  were  sampled 
at  an  interval  of  1.2  msec.  Since  we  are  mainly  interested  in  the 
alpha  component  of  the  EEG,  the  data  were  filtered  to  obtain  the 
components  around  the  alpha  range  (7-13  Hz)  through  a transversal  filter. 
The  computer  program  reads  this  EEG  data  as  its  input.  Other  required 
information  such  as  the  number  of  sample  points  N for  which  the  process 
is  considered  stationary,  the  order  of  the  process  P , the  resampling 
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value  K , the  pr< diction  lead  length  L and  the  window  width  NWND 
for  updating  the  parameters  are  also  read  in  as  input.  (For  efficiency, 
the  parameters  are  not  updated  with  each  new  data  point,  but  only  after 
NWND  new  data  points  have  been  read  in.)  A flow  chart  of  the  program 
is  shown  in  Figure  33.  As  soon  as  the  above  variables  and  the  initial 
N points  are  entered,  the  program  computes  the  least-square  estimates 

of  the  parameters  by  solving  the  set  of  linear  equations  (7).  The  error 

, 2 

variance  ag  is  also  estimated.  These  estimates  are  then  used  to 
generate  minimum  mean-square  error  forecasts  up  to  L steps  ahead; 
i.e.,  forecast  of  t + KL  is  obtained  at  time  t . The  variance  of 
the  forecast  errors  and  the  confidence  interval  for  a certain  significance 
level  are  also  calculated.  At  this  time,  the  next  data  point  is  read 
in;  the  residual  is  calculated  as  the  difference  between  the  actual 
value  and  the  forecast;  and  the  forecasts  for  the  remaining  points  are 
updated  by  using  Eq.  (5).  A new  forecast  for  the  value  L steps 
ahead  is  produced.  This  process  is  repeated  until  NWND  new  points 
have  been  exhausted.  The  procedure  will  repeat  again  starting  with 
the  recalculation  of  the  parameter  estimates. 

Results  of  the  prediction  for  various  lead  times  are  shown  in 
i igures  34  through  38.  In  order  to  compare  the  predicted  waveforms  with 
the  actual  EEG,  we  plotted  the  predicted  waveforms  and  then  superimposed 
the  actual  EEG  on  top  of  them.  The  order  of  the  autoregressive  process 
was  chosen  as  p = 7 by  examining  the  partial  autocorrelation  function 
of  EEG  data  for  k = 10  . The  effective  lead  time  in  seconds  is 
kiT  where  T is  the  sampling  interval  and  i denotes  the  lead  time 
in  number  of  samples.  In  Figures  34  through  38,  the  waveforms  were 
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Figure  33  A flowchart  for  the  realization 
of  the  prediction  scheme 


figure  34  Actual  and  predicted  waveform 
for  lead  time  of  24  msec. 


Figure  35  Actual  and  predicted  waveform 
for  lead  time  of  48  msec. 


w 
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igure  39  Prediction  error  variance 
vs.  lead  time 
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plotted  against  the  sample  numbers  and  the  lead  time  i was  indicated 
at  the  lower  right  corner.  For  instance,  in  Figure  36,  £ = 5 means 
a lead  time  of  10  X 5 X 1.2  = 60  msec.,  sin'-e  k = 10  and  T = 1.2  msec; 
hence,  the  length  indicated  for  £ = 5 represents  60  msec.  The  amplitude 
scale  is  relative.  The  solid  curves  represent  the  actual  EEG  and  the 
dotted  curves  are  the  predicted  waveforms.  In  Figure  34,  we  see  that 
for  a lead  time  of  24  msec.,  the  two  curves  practically  coincide  with 
each  other.  This  indicates  a prediction  with  little  error.  We  progressive- 
ly increased  the  lead  time  from  48  msec,  in  Figure  35  to  120  msec, 
in  Figure  38.  It  is  observed  from  this  series  of  figures  that  the 
coincidence  of  the  predicted  and  the  actual  waveforms  worsens.  This 
indication  of  increasing  error  as  lead  time  increases  is  expected, 
lor  a quantitative  comparison,  we  computed  the  prediction  error  variance, 
which  is  used  as  a performance  measure  of  the  prediction  scheme,  for 

various  lead  times.  This  forecast  error  variance  is  shown  in  Figure  39. 

2 

The  vertical  scale  gives  VU)/E(xt>  which  represents  the  prediction 
error  variance  relative  to  the  variance  of  signal  to  be  predicted, 
lor  example,  for  £ = 5 (i.e.,  a lead  time  of  60  msec.)  the  error 
variance  is  1$  of  the  signal  variance.  This  gives  a measure  of  how 
close  the  predicted  values  will  be  to  the  actual  values.  We  can  see 
at  a glance  from  Figure  39  that  the  performance  deteriorates  as  let^d 
time  increases.  In  using  this  scheme,  one  can  determine  a maximum 
lead  time  for  any  given  tolerance  level.  For  our  purposes,  this 
scheme  is  satisfactory. 
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IV.  CHAUACTFHIZATION  OF  SCAN  PATTERNS 


One  ol  the  important  aspects  of  the  research  done  in  this  project 
has  been  the  development  of  techniques  and  tools  for  the  automatic 
analysis,  characterization  and  ri^play  of  scan  patterns.  In  visual 
tasks,  the  sequence  of  fixations  usee’  in  scanning  the  visual  target 
is  called  the  "scanpath",  and  we  will  refer  to  a "scan  pattern"  as  a 
collection  of  information  about  the  properties  of  one  or  more  scan- 
paths  recorded  by  one  or  several  different  observers.  Naturally, 
to  be  able  to  differentiate  a superior  scan  pattern  which  results  in 
superior  visual  memory  from  an  inferior  one,  we  must  first  be  able  to 
characterize  a scanpath  in  some  terms  other  than  just  a list  of  (X,Y) 
fixation  coordinates, 

A graphical  illustration  of  a sequence  of  the  foveal  fixations 
is  shown  in  figures  40  and  41.  figure  40  shows  the  outline  drawing 
of  a still  life  scene  which  was  inspected  by  a subject  with  his  scanpath 
superimposed.  In  figure  41  we  have  plotted  a sketch  ol  the  visual 
stimulus  impressed  on  the  fovea  for  each  fixation  point  shown  in 
figure  40.  The  sketches  show  a visual  field  of  6°  in  diameter,  a 
region  which  includes  all  of  the  fovea  under  practically  any  definition. 
I hese  drawings,  then,  represent  a time  sequence  ol  snapshots  which  are 
the  visual  input  to  the  observer  who  must  create  a unified  picture  of 
the  scene  being  viewed  1 rom  the  sequence  of  picture  fragments  shown  in 
the  sketches. 

One  observation  which  has  often  been  made  is  that  the  fixation 
points  tend  to  be  grouped  around  certain  features  or  areas  of  the 
visual  image,  and  are  relatively  sparse  elsewhere.  Thus,  we  can  measure 
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Figure  40  Computer  derived  scanpath  superimposed 
on  sketch  of  original  scene 
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Figure  41 


Sequence  of  foveal  fixations 


the  nearness  of  one  fixation  point  to  another  and  assign  a similarity 
rioasure  to  the  two  points  based  upon  this  distance.  The  ability  to 
group  or  cluster  a set  of  fixation  points  is  an  important  step  in  the 
characterization  of  a scanpath,  and  several  different  techniques  have 
been  used  to  measure  point  similarities  and  to  assign  points  to  clusters. 

The  use  of  a "point  similarity  function"  for  clustering  fixation 
points  has  been  reported  earlier  and  will  not  be  reviewed  here. 

We  have  also  used  two  alternate  algorithms  for  selecting  cluster 
centers,  each  with  its  own  particular  advantages  and  disadvantages. 

These  algorithms  are  known  as  "ISODATA"  and  "Minimum-Spanning-Tree-Clustering" 
or  MSTC. 

In  the  ISODATA  method,  the  number  of  clusters  desired  is  given 

an  input  variable,  and  the  algorithm  partitions  the  fixation  points 
into  subsets  such  that  the  total  distance  between  the  fixation  points 

in  a cluster  and  the  cluster  center  is  minimized  for  all  clusters. 

The  cluster  center  starting  points  are  usually  selected  randomly,  and 
the  algorithm  keeps  moving  the  cluster  centers  until  the  above 
criterion  has  been  reached. 

The  MSTC  algorithm  automatically  determine  the  number  of  clusters 
in  the  data  and  assigns  fixation  points  to  the  clusters.  Three 
variables  given  as  input  may  be  adjusted  to  change  the  criteria  by 
which  clusters  are  selected. 

The  chief  advantages  or  disadvantages  of  these  two  algorithms  are 
summarized  below: 
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I SODATA 


(1)  Required  the  desired  number  of  clusters,  which  usually  is  not 
known. 

(2)  Solutions  often  are  not  unique,  but  vary  with  starting  conditions. 

(3)  The  iterations  can  require  a substantial  amount  of  computing. 

(4)  Many  points  can  be  easily  accommodated. 

MSTC 

(1)  The  solution  is  unique,  and  the  number  of  clusters  is  automatically 
determined. 

(2)  The  compu  ation  is  quite  fast. 

(3)  Many  points  requires  a large  amount  of  storage. 

Our  approach  to  characterising  scan  patterns  is  based  on 
the  following  two  assumptions: 

(a)  There  are  a discrete  number  of  "lixation  centers"  in  the  scene 
being  viewed  and  each  fixation  point  in  the  scan  can  be  assigned 

to  one  oi  the  fixation  centers.  This  process  of  assignment  or 
clustering  is  done  using  the  MSTC  algorithm  described  above.  The 
center  of  the  clusters  should  be  in  close  agreement  with  the 
centers  of  the  fixation.  Due  to  the  random  nature  of  the  fixation 
points,  in  any  repetition  (cycle)  of  the  scanpath,  the  fixation 
points  do  not  always  coincide. 

(b)  A scanpath  then  becomes  a sequence  of  transitions  from  one 
cluster  center  to  another.  Any  saccade  with  starting  and  ending 
points  in  the  same  cluster  is  discarded. 
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T rom  the  set  of  all  saccades  leaving  a given  cluster  center,  we 
can  find  the  transition  or  transitions  which  have  the  highest  probability. 
These  transtions  are  termed  the  "most-probable-saccades".  ; he  "most- 
probable-scanpath"  is  then  defined  to  be  the  set  of  all  most-probable- 
saccades.  The  most-probable-scanpath  depends  on  how  the  clusters  arc 
chosen,  and  is  not  necessarily  a closed  path  through  the  cluster  centers. 
Its  usefulness  is  that  it  concisely  summarizes  the  most  important  aspects 
of  the  scanpath  used  by  the  observer  in  viewing  the  given  visual  target. 

A scan  pattern  is  some  set  of  most-probable-scanpaths  for  a given  observer 
viewing  the  same  or  different  targets,  or  for  several  observers  viewing 
the  same  target,  etc. 

figures  42  and  43  illustrate  how  these  methods  may  be  employed  to 
reduce  scanpath  data  to  a more  concise  form.  Figure  42  is  a graph  of  a 
scanpath  superimposed  upon  a line  drawing  representing  the  still  life 
scene  which  the  observer  was  viewing.  The  first  recorded  fixation 
(1)  is  in  the  center  of  the  avocado.  The  last  fixation  (45)  is  in  the 
middle  of  the  carrot.  The  viewing  time  required  to  create  this 
scan  totaled  20  seconds.  Figure  43  represents  the  same  data  after 
processing  by  the  MSTC  program.  The  program  partitioned  the  fixation 
points  into  11  subsets  or  clusters  (using  certain  criteria  supplied  by 
the  programmer)  and  then  computed  the  most  likely  transitions  from 
cluster  to  cluster.  Clusters  containing  only  one  fixation  point  were 
ignored.  The  resulting  plot  shows  the  most  likely  saccades  between  areas 
containing  large  numbers  of  fixations,  and  presents  the  original  data 
in  a highly  condensed  format. 

In  order  to  compare  one  set  of  fixations  or  one  scanpattern 
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;ure  42  Scan  of  still  life  scene 
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against  another,  criteria  must  be  established  which  can  be  used  to 
evaluate  subtle  differences  in  the  following  characteristics: 

(1)  Total  numuer  of  clusters  - The  MSTC  algorithm  partitions 
all  the  lixation  points  into  some  number  of  clusters,  very 
rarely  will  the  number  of  clusters  in  two  different  scans  be 
the  same.  Rules  must  be  established  for  eliminating  or  merging 
clusters  which  complicate  the  analysis.  for  example,  clusters 
containing  a single  point  can  usually  be  removed. 

(2)  Similarity  between  clusters  - The  geometrical  center  of 
two  clusters  will  . _.ually  be  different,  and  yet  it  may  seem 
obvious  that  both  clusters  are  associated  with  or  centered  on 
some  feature  of  the  picture. 

(S)  The  most  probable  saccades  between  clusters  will  often  be 
different  in  two  scans,  even  if  the  cluster  centers  are  nearly 
alike.  The  similarity  ol  two  scans  is  judged  by  their  transition 
probabilities. 


The  above  criteria  are  dependent  upon  the  type  of  experiment  being  run 
and  no  1 ixed  set  of  rules  will  be  suggested  here.  Based  on  the  above 
criteria  and  using  statistical  methods,  we  are  able  to  show  from  our 
preliminary  results  that  the  scan  patterns  of  the  same  visual  material 
by  the  same  subject  who  scored  very  high  on  Marks'  Visual-Memory  Task 
are  consistent. 
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The  measurement  and  prediction  oi  eye  movements  and  related  1.10 
potentials  are  necessary  prerequisites  lor  the  control  of  visually 
displayed  material  which  can  enhance  an  observer’s  perception  or  memory 
or  a scene  or  visual  environment.  This  s it  ion  describes  some  means 
to  couple  the  measurement  oi  eye  movements  with  the  control  ol  a 

computer-driven  graphical  display. 

The  task  undertaken  was  to  provide  a lacility  lor  dynamically 
controlling  the  loveal  or  peripheral  field  of  view  ol  an  observer.  Such 
a lacility  would  provide  a flexible  method  lor  studying  interaction 
between  eye  position  and  changes  in  the  visual  field.  Some  of  the  ob- 
jective’s ol  this  research  are; 

(a)  to  develop  techniques  lor  dynamically  controlling  the  visual 
field  in  humans  without  any  restrictive  mechanical  attachment  to 

the  eye; 

(b)  to  obtain  a subjective  feeling  for  perception  without  use  ol 
the  noimal  visual  field;  and 

(c)  to  quantitatively  measure  changes  in  the  observer’s  performance 

(and  scanpath  characteristics)  as  a function  of  field  size  and  type. 
The  system  works  as  follows:  At  the  start  ol  each  refresh  of  the 


display  (30  times  per  second),  the  position  of  the  eye  is  measured  and 
the  corresponding  point  on  the  display  is  computed.  Parameter,  within 
the  program  determine  v;hat  portion  ol  the  visual  field  will  be  displayed 
lovea  only,  periphery  only,  and  the  size  of  the  field  in  either  ease. 
Then,  each  point  in  the  display  list  is  checked,  and  il  it  is  within  the 
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required  Held,  it  is  displayed.  Figures  44  and  45  depict  the  stimulus 
for  the  two  modes  of  display.  The  box  (which  was  used  instead  ol  a 
circle  to  simplify  computations)  remains  centered  on  the  loven,  no  matter 
where  on  the  screen  the  observer  looks. 

The  visual  scenes  used  in  these  experiments  were  outline  drawings 
of  simple,  still  scenes.  No  text  or  highly  detailed  pictures  were 
involved.  A new  scene  could  be  entered  into  the  computer  through  use 
of  a graf/pen  digitizer.  Once  entered,  the  points  in  the  picture  wet e 
stored  by  X,Y  position,  and  a display  list  was  created  for  later  use 
by  the  real-time  interactive  program. 

Alter  adjusting  the  eye  tracker  and  calibrating  the  instrument, 
the  observer  would  be  presented  with  a field  of  view  which  could  be 
either  a fovea-only  presentation  or  a periphery-only  presentation. 

The  size  of  the  field  (width  of  box  boundary)  was  also  adjustable,  i rom 
less  than  one  degree  to  the  maximum  width  of  the  display.  Use  of  Lhi*- 
program  has  produced,  so  far,  these  results: 

Impaired  perception  for  fovea-only  vision  - If  only  a two  or  three  degree 
wide  foveal  field  is  displayed  (when  the  display  itself  spans  twenty 
or  more  degrees),  the  observer  finds  it  difficult,  if  not  impossible,  to 
perceive  the  display.  What  the  observer  sees  are  short  line  segments, 
and  sometimes  junctions  between  line  "egnmuts.  The  only  way  to  perceive 
an  object  is  to  consciously  follow  the  contours  and  mentally  reconstruct 
the  path  followed.  Simple  objects  take  several  minutes  to  detect.  Few 
observers  get  a complete  idea  of  the  entire  scene.  Far  more  use  of 
peripheral  visual  information  is  made  than  was  anticipated.  Large 
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saccades  must  be  forced,  and  the  observer  becomes  very  conscious  of  having 
to  manipulate  his  own  eye  movements. 

Perception  with  peripheral  vision  - If  the  foveal  region  is  blanked, 
little  perceptual  loss  occurs  until  the  window  size  approaches  roughly 
10  , that  is,  the  scene  is  still  perceived  even  when  the  window  is 
Quite  large.  Our  stimuli  did  not  include  text  or  other  material  where 
high  visual  acuity  was  essential.  This  result  is  somewhat  .«  uprising 
in  view  of  the  general  assumption  of  the  importance  of  foveal-hign 
acuity  field.  The  peripheral  field  apparently  makes  an  important 
contribution  to  the  overall  perception  of  a scene,  even  the  visual 
acuity  may  be  quite  low  in  the  outer  regions. 

Importance  of  closed  contours  - If  the  fovea-only  visual  field  is  being 
displayed,  little  improvement  in  perception  of  objects  within  the 
field  occurs  until  the  window  size  is  large  enough  to  include  the  complete, 
closed  contour  of  an  object.  The  absolute  size  of  the  foveal  region  seems 
to  be  less  importrnt  than  the  fact  that  complete  objects  are  visible, 
at  least  for  the  line  drawings  we  used.  This  result  supports  another 
observation  we  have  made,  which  is  that  fixations  tend  to  be  centered 
on  the  main  surfaces  of  a simple  object  and  are  not  as  often  bound  on 
the  edges  or  corners  of  the  object. 

Appearance  of  a phantom  window  - After  about  ten  seconds  of  viewing,  and 
especially  for  the  fovea-only  condition,  a phantom  window  with  pronounced 
grey-black  edges  appears.  The  subject  sees  a grey  background  with  a 
black  window  cut  in  the  center,  through  which  the  white  lines  may  be 
viewed.  The  window,  of  course,  is  fixated  on  his  eye  like  an  afterimage. 
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Two  causes  for  the  window  are  evident.  Only  the  foveal  region  receives 
stimulation,  and  the  rods  and  cones  adapt  to  a higher  ambient  light 
level  than  in  the  periphery,  which  adapts  to  "visual  grey".  Also,  the 
sharp  truncation  of  many  contours  along  a straight  line  reinforces  the 
perception  of  an  edge.  Similar  effects  have  been  noted  before  by  Yarbus. 

These  preliminary  experiments  with  an  eye-movement  controlled 
visual  display  have  proved  quite  interesting  and  useful.  The  next  step 
is  to  incorporate  the  predictive  algorithms  discussed  in  other  previous 
sections  of  this  report  to  present  the  visual  cues  at  the  desired 
locations  of  the  visual  field  at  the  desired  moments. 
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VI.  CONCLUSIONS 


Our  aim  has  been  to  develop  and  use  biocybernetic  techniques  for  the 
enhancement  of  visual  memory.  To  this  end,  we  have  designed  a closely- 
coupled  man-machine  system  f-<r  implementation.  We  have  developed  and 
completed  the  models  for  saccadic  eye  movements  and  KKC  signals  in  this 
man-machine  system.  The  models  are  used  for  monitoring  and  prediction  of 
eye  positions  and  brain  states.  With  the  continuously  up-dated  information 
concerning  the  eye  position  and  brain  state  for  adjusting  the  stimulus 
parameters  and  the  monitoring  and  prediction  schemes,  we  can  guide  the 
eyes  to  fixate  at  the  specified  locations  of  the  visual  field  at  the 
specified  instants  of  time  through  a technique  we  developed.  If  these 
specified  locations  and  instants  are  the  optimally  required  for  superior 
visual  memory,  then  the  vividness  and  persistence  of  the  desired 
after-image  will  be  enhanced.  We  have  approached  the  problem  of 
determination  of  these  optimal  locations  and  time  instants  by  monitoring 
eye  movements  and  LEG  signals  through  the  models  which  we  have  developed. 

Through  the  EEG  model,  we  have  obtained  a better  assessment  of  EEC's 
role  as  the  timing  mechanism  for  visual  information  acquisition  and 
processing.  We  have  also  shown  that  the  scan  pattern  is  more  consistent  for 
individuals  with  good  visual  men ory  (as  scored  by  Marks'  Visual-Memory 
Task)  than  those  with  poor  visual  memory.  This  comparison  is  made  possible 
by  our  method  of  characterization  of  scan  patterns.  With  these  results, 
it  is  possible  to  determine  the  optimal  locations  in  the  visual  field 
and  the  optimal  time  instants  for  presenting  visual  stimuli.  We  conclude 
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that  the  approach  for  visual  memory  enhancement  as  proposed  is  a feasible 
one.  At  the  termination  of  this  project,  we  have  implemented  all  of 
the  parts  shown  in  the  closely-coupled  man-machine  system  but  the 
feedback  path. 

Significant  accomplishments  made  during  this  research  have  been 
published  or  presented  in  international  or  national  scientific  conferences 
and  two  Ph.D.  dissertations  are  near  completion. 
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appendix  a 


ELECTRONIC  INTERFACE  SYSTEM 


**  ShaM  de”rlbe  here  ,he  interface  system  for  eye-movement  data 
acquisition  only.  We  ere  pre.ently  using  a double  Purkinje  image  eye 
tracker  developed  a,  Stanford  Research  Institute  (SRI)  for  measuring 
eye  movements.  I,  relies  on  measuring  the  motion  of  the  reflections 
from  the  front  surface  of  the  cornea  and  from  the  hack  surface  of  the 
lens  of  the  eye  (the  first  and  fourth  Purkinje  image  ,.  The  instrument 
is  discussed  in  an  article:  "accurate  two-dimensional  eye  tracker 
using  first  and  fourth  Purkinje  images",  by  T.  S.  Cornsweet  and  H.  D. 
Cra”C  *he  i°ynal  Of  the  Optical  Society  of  Amer.e,  vol.  63_  8> 

PP.  921-928,  August,  1973.  figures  A.l  and  a. 2 show  the  set  up  of  this 
eye  tracker. 

The  first  and  fourth  Purkinje  images  are  generated  by  positioning 
a narrow  beam  of  infrared  light  on  the  subject’s  pupil.  The  image  of 
the  eye  and  its  attendant  Purkinje  image  reflect  off  the  infrared  mirror 
in  front  of  the  subject,  through  two  large  collimating  lenses  and  on  to 
a movable  mirror,  a four-quadrant  photodetector  senses  the  position  of 
the  first  Purkinje  image.  Signals  from  these  four  quadrants  are  used 
to  drive  two  high-speed  servo  motors  to  detect  the  movable  mirror  in 
altitude  and  aalmuth.  The  function  of  this  servo  system  is  to  position 
the  image  of  corneal  reflection  when  the  eye  moves  so  that  the  reflection 
is  always  in  the  same  position;  thus,  a stable  reference  of  the  eye  is 
provided,  as  the  eye  moves,  the  image  reflected  by  the  movable  mirror 
•ill  remain  stationary,  a second  optical  system  in  tandem  consisting  of 
another  movable  mirror  and  quadrant  photodetector  tracks  the  fourth 
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Figure  A.l  View  of  subject  in  the  head 
positioner  with  eye-tracker  electronics 


Figure  A. 2 View  of  eye-tracker 
optics  and  subject 


ll 

1 ■ 

1 D ' ■' 

B* 

jg 


Purkinje  image  and  measures  its  movement  relative  to  the  first  image. 

The  two-dimensional  motion  of  this  image  derived  from  the  position  of 
this  mirror  system  is  used  as  the  output  signal  of  the  subject’s  eye 
movements,  The  system  is  designed  to  measure  eye  movements  with  an 
accuracy  of  up  to  two  minutes  of  arc.  The  overall  system  has  a flat 
frequency  response  of  up  to  100  Hz. 

To  assure  the  proper  use  of  the  eye  tracker  requires  elaborate 
procedures  for  alignment  of  subjects  and  calibration  of  the  eye  tracker. 
The  calibration  is  done  through  a computer  program  (MAP)  we  developed. 
This  calibration  program  serves  to: 

(1)  provide  a way  of  relating  the  output  voltages  of  the  instrument 
to  the  stimulus  position*; 

(2)  give  an  estimate  of  the  size  of  the  visual  field  over  which 
the  instrument  is  able  to  irac  ; the  eye  movements;  and 

(3)  indicate  whether  the  eye  tracker  is  correctly  adjusted. 

It  has  been  a major  problem  to  relate  the  measured  eye  positions 

to  the  actual  locations  on  the  scene.  This  difficulty  arises  mainly 
from  the  distortion  introduced  in  measuring  eye  movements.  To  alleviate 
this  difficulty,  we  have  developed  and  implemented  a scheme  which  maps 
automatically  the  eye  fixations  onto  the  scene  without  distortion.  This 
is  achieved  by  expressing  both  the  fixations  and  the  line  drawing  of 
the  inspected  scene  in  terms  of  visual  degrees  relative  to  a calibration 
slide.  Since  both  are  on  the  same  coordinate  system,  they  can  be  easily 
superimposed  without  any  fear  of  distertion. 

Our  scheme  consists  of  the  following  computer  programs:  MAP, 

HXFIL,  OVLAY,  and  PFIX.  MAP  produces  a calibration  file  relating 
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eye  tracker  voltage  with  the  visual  angle  subtended  by  the  eye.  FIX  ML 
analyzes  tne  data  collected  by  the  DIGIT  program.  It  reduces  the  data 
into  fixations  and,  with  the  aid  of  the  calibration  file  produced  by 
MAP,  it  translates  the  location  of  these  fixations  into  units  of  visual 
angle.  The  origin  of  this  visual  coordinate  system  is  assumed  to  be 
the  center  point  of  the  calibration  slide.  OVLAY  is  a variation  of  the 
DRAW  program  which  uses  the  graf/pen  to  digitize  pictures.  OVLAY  is 
used  to  digitize  an  outline  of  the  inspected  scene  on  a slide.  Some 
additional  data  are  added  to  the  file  to  enable  it  to  be  eventually 
translated  into  units  of  visual  angle  with  the  exact  coordinate  system 
by  HXHL.  PI  IX  plots  the  fixations  and  superimpose  the  outline  of 
the  scene.  This  can  easily  be  done  since  the  fixation  locations  and  the 
outline  are  expressed  in  the  same  visual  angle  coordinate  system. 

To  illustrate  the  effectiveness  of  this  scheme,  we  show  an  example 
in  figures  A. 3 and  A. 4.  The  observer  was  asked  to  fixate  at  those 
circled  points  of  a scene  as  shown  in  Figure  A. 3.  The  eye-movement  data 
taken  were  then  processed  by  the  mapping  scheme.  The  result  was  plotted 
as  shown  in  Figure  A .4  . it  is  seen  that  the  measured  eye  fixations 

correspond  very  closely  co  those  circled  points  in  Figure  A. 2 in  the 
right  sequence. 

It  should  be  emphasized  that  the  process  described  is  fully  automatic 
and  will  compensate  for  different  subjects  r.nd  any  of  the  permissible 
variations  in  the  experimental  set  up  such  as  size  of  stimulus,  distance 
between  the  screen  and  the  observer,  etc.  With  this  system,  we  are 

capable  of  processing  massive  amounts  of  eye-movement  data  with  relative 
ease. 
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Figure  A. 3 Requested  fixation  points 
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Figure  A. 4 Responses  to  req jested  fixation  points 
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